Vanishing neoclassical viscosity and physics of the shear layer in stellarators 
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The drift kinetic equation is solved for low density TJ-II plasmas employing slowly varying, time- 
dependent profiles. This allows to simulate density ramp-up experiments and describe from first 
principles the formation and physics of the radial electric field shear layer. The main features of the 
transition are perfectly captured by the calculation, and good quantitative agreement is also found. 
The results presented here, that should be valid for other non-quasisymmetric stellarators, provide 
a fundamental explanation for a wealth of experimental observations connected to the shear layer 
emergence in TJ-II. The key quantity is the neoclassical viscosity, which is shown to go smoothly 
to zero when the critical density is approached from below. This makes it possible for turbulence- 
related phenomena, and particularly zonal flows, to arise in the neighborhood of the transition. 



Understanding transport barriers in magnetic confine- 
ment devices [1] is a crucial issue for the fusion pro- 
gramme because they allow access to regimes of improved 
confinement in which commercial fusion reactors may be 
viable. Confinement transitions have been observed and 
comprehensively documented in a large variety of toka- 
maks and stellarators (see Ref. [2] and references therein). 
While the reduction of turbulence by sheared flows is gen- 
erally accepted as a key ingredient in the transitions [3] , 
the identification of the physical mechanisms that create 
these sheared flows is still an open issue. Zonal flows 
(i.e. turbulence-generated, fluctuating and flux-surface 
collective radial electric fields) are considered a possi- 
ble catalyst of the transition, since they are observed 
close before confinement transitions and they are able to 
regulate transport [4]. In non-quasisymmetric stellara- 
tors, the non-ambipolar radial fluxes, and consequently 
the radial electric field, are basically determined by the 
neoclassical theory of collisional transport in magnetized 
plasma [5]. Therefore, the interaction between neoclassi- 
cal and turbulent processes, and concretely the extent to 
which turbulence can overcome the neoclassical viscos- 
ity and modify the E x B rotation through momentum 
transport, has been the subject of recent works [4-7]. 
In this Letter, we contribute to the above programme 
with the study of the so-called low density confinement 
transition [8-12] in the flexible heliac TJ-II [13] which, 
as we discuss below, shows similarities with other tran- 
sitions described in the literature. In this work the for- 
mation of the shear layer is described from first prin- 
ciples in the framework of neoclassical theory and the 
turbulent phenomena that arise in the neighborhood of 
the transition are shown to be regulated by neoclassi- 
cal transport. TJ-II [13] undergoes a spontaneous con- 
finement transition typically at a line-averaged electron 
density n^ = n cr ~0.6 x 10 19 m -3 . At this empirical crit- 
ical density n cr , the radial electric field E r changes from 
positive to negative and, at the same time, a transport 
barrier is generated close to the edge. The reversal of E r 
starts where the density gradient is maximum (p « 0.8, 
where p = r/a is the normalized radius) and then prop- 



agates across the entire region 0.5 < p < 0.9 at a speed 
of the order of several m/s [14]. A sheared E r always 
appears around the ^-reversal point. Further increase 
in does not modify qualitatively E r . It is argued [10] 
that the transport barrier is caused by turbulence reduc- 
tion when the shearing rate of the E x B flow reaches 
values of the order of 10 5 s _1 ; that is, when it becomes 
comparable to the estimated linear growth rate of the 
resistive interchange instability in TJ-II [15]. 

Whereas similar transport barriers, related to jumps 
between roots of the ambipolar equation, have been ob- 
served in the low density regime of other stellarators [16], 
the very detailed study carried out in TJ-II has re- 
vealed additional interesting phenomena during the tran- 
sition [9, 11]. First of all, the level of turbulence and 
the E x B flux are seen to increase prior to the transi- 
tion, and this was later shown to be associated to long- 
range-correlated (LRCed) electrostatic potential struc- 
tures that grow when approaching the critical density. 
These experimental results have historically led to an in- 
terpretation of the transition in terms of the paradigm 
of sheared flow generation by Reynolds stress and tur- 
bulence reduction by shear, see e.g. Refs. [8, 10]. Along 
the same line of reasoning, the peaking of the potential 
relaxation time in biasing turn-off experiments [10] and 
of the shear-flow susceptibility in electrode-biasing ex- 
periments [17] were interpreted in terms of an increased 
turbulent drive during the transition. The increase of 
the level of turbulence and of the long-range correla- 
tions are neither specific of this transition nor of stel- 
larators: they are also observed in the L-H transition of 
TJ-II (which happens at higher n£) and of other devices, 
see e.g. Ref. [18] and references therein. We note that 
there are important differences between the low-density 
and the L-H transition (e.g. the improvement of energy 
confinement in the latter). Nevertheless, the neoclassical 
modelling of the L-H transition [19] also relies on bifur- 
cations of the ambipolarity condition. 

Since the inversion of E r is an essential part of the 
transition, a neoclassical study of the transition is in or- 
der. The existence of multiple roots of the ambipolar 
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equation is a well-known feature of stellarators in the 
low-density regime [20]. Indeed, a static neoclassical cal- 
culation of E r for plasmas immediately below (above) n cr 
was performed in Ref. [21], yielding positive (negative) 
E r in qualitative agreement with the experiment. In this 
work, we go beyond previous calculations and perform a 
full dynamical neoclassical calculation of the formation 
of the sheared E r . We simulate the density ramp- up 
that leads to the transition and describe, from first prin- 
ciples, the formation and evolution of the shear flow in 
good agreement with the experiment. Furthermore, we 
show for the first time that the behaviour of the three 
quantities discussed above (amplitude of low-frequency 
LRCed potential fluctuations, potential relaxation time 
and shear-flow susceptibility) appear here as natural con- 
sequences of a neoclassical bifurcation. Namely, as the 
change of root approaches, the neoclassical poloidal vis- 
cosity (the restoring force of E r deviations towards its 
ambipolar value), goes to zero. It will be shown that this 
automatically implies a maximum of the low-frequency 
plasma potential fluctuations, of the relaxation times and 
of the shear-flow susceptibility. 

Let us begin by deriving the equation for the electric 
field evolution. We start from the momentum balance 
equation [5] summed over species: 

m i ^^+V-n i +V-n e =jxB. (1) 

Here, u is the ion flow tangent to flux surfaces, II S is the 
viscosity tensor, U s = m s j vv/ s (x, v, t)d 3 v, or momen- 
tum flux of species s, f s its distribution function, and 
j x B is the Lorentz force. Other forces (such as neutral 
friction) may be included as extra terms on the RHS of 
Eq. (1). We have assumed a quasineutral plasma con- 
sisting of singly charged ions and electrons (n e = n^ = n). 
Note that we have dropped the inertia of the electrons, 
given their much lower mass m s , but kept the electron 
viscosity tensor as it cannot be neglected in our low-n, 
high-T e plasmas [22]. We work in Hamada magnetic co- 
ordinates (-0,0, £), and follow the notation of Ref. [6]. 
The lowest order incompressible ion flow is conveniently 
written as: 

u = 2 7T (^±+<j ) '(iP))e +A(i;)B. (2) 
\ ne J 

The flux surface label ip is the toroidal magnetic flux 
and the prime stands for derivative. The first term 
on the RHS of Eq. (2) contains the diamagnetic and 
E x B perpendicular flows (pi is the ion pressure, e the 
elementary charge, <p is the electrostatic potential, and 
e<9 x B = (27r) _1 VV ; ) together with the parallel Pfirsch- 
Schliiter flow (V • e# = and (e# • B) = for a cur- 
rentless stellarator) . The term AB is the ion bootstrap 
flow [23]. If we project Eq. (1) along eg and take flux- 
surface- aver age (•) we obtain our evolution equation for 



the radial electric field: 

dEr = ld_ (ri(r)\ _ E ldn + 
dt n dt \ e J r n dt 

+ A -r (e(r e -rO + (j-Vr)), (3) 

47r z mn{ee • ee) 

where E r = — <//(r) and the minor radius r is a geo- 
metric flux label defined in terms of the volume V(r) = 
7rr 2 L ax , where L ax is the length of the magnetic axis. 
We have obtained the radial particle fluxes from T s = 
— - |tt^ (eg • V • U s ) [24]. The viscosity tensor can be 
split into a neoclassical part, given by the gyrotropic pres- 
sure tensor, and an anomalous contribution: 

n s = nf c + nr = p s \\hb + Ps± (i - bb) + u a s n . (4) 

As mentioned above, in non-quasisymmetric confining 
magnetic topologies [5], the leading order contribution 
to Eq. (3) is (e<9 • V • n^ c ), being much larger than 
(eg • V • n^ n ), which will be therefore neglected. Higher 
order neoclassical terms like the shear-flow viscosity [5] 
can be considered. This term arises from orbit deviations 
away from flux surfaces and transforms Eq. (3) into a 
non-local diffusion equation for E r . Since this compli- 
cates the discussion and does not fundamentally modify 
the predictions [26], we neglect it and provide a posteriori 
justification based on the experimental results. 

We use the Drift Kinetic Equation Solver (DKES), 
complemented with momentum-correction techniques, to 
evaluate the pressure anisotropy in the TJ-II magnetic 
field in the parameter range usually found experimen- 
tally in the vicinity of the transition. Details of the 
calculation and convolution of the monoenergetic coef- 
ficients may be found in Ref. [23] and references therein. 
From Eq. (3), the time-evolution of E r is fully deter- 
mined if we know, at every instant of time, the magnetic 
configuration b mn and the profiles n, T e , T^. Since we 
simulate a pure proton-electron plasma [22], the effec- 
tive charge Z e ff (which mainly affects collisionality) is 
set equal to one. We perform a numerical simulation of 
a density ramp across the critical density for a plasma 
with profiles n(p,t) : T e (p,i) and Ti(p,t) that mimic the 
experimental ones, see Fig. 1 and Refs. [9, 21]. We set 
(j • Vr) = unless otherwise stated. This is implied by 
quasineutrality (V • j = 0), but a net radial plasma cur- 
rent can be induced in plasma biasing experiments. An 
important point to be noted is that, although we assume 
that the leading non-ambipolar particle fluxes are neo- 
classical, we make no particular assumption on the total 
particle or energy fluxes. They are included (together 
with the sources) implicitly in the evolution of n, T e , T^. 

In Fig. 2 we sketch the evolution of the main quantities 
of our numerical ramp-up experiment in a representative 
radial position, p = 0.7. The evolution of E r is given 
only by the evolution of n and T e during the ramp-up. 
As n rises, the radial fluxes become larger and E r less 
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FIG. 1: (Color online) Plasma radial profiles for 
selected times: low (high) n in open (closed) circles. 
The evolution is given through ^§r=3s _1 , 



^^ = -1.5s- 1 ,and^=0. 
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FIG. 2: (Color online) Time evolution of the relevant 
quantities at p = 0.7 during the density ramp-up. 



positive. is slightly higher than T e , which sets the sign 
of dE r /dt, but the difference is not even visible in the 
figure. This is a consequence of the very slow variation 
of n and T e , and will allow us to discuss some of the 
results in terms of the steady-state ambipolar equation. 
For n~0.6 x 10 19 m -3 , there is a change of root: E r goes 
from positive to negative in several tens of ps. For larger 
n, further increase in n leads to smaller Y s and more 
negative E r . This change of behaviour of T s (they grow 
with n below n cr and decrease above n cr ) might indicate 
a neoclassical particle transport barrier (nevertheless, let 
us remember that it is probably the ambipolar anomalous 
flux [9] which controls particle transport). 

This general behaviour is expected for neoclassical sim- 
ulations of TJ-II low-density plasmas [21, 22]. In Fig. 3, 
we show the ambipolar equation at p = 0.7 for several 
relevant times. Since we start from low collisionality, E r 
is positive (t = 10 ms). As n is raised, a negative stable 
root appears (t = 50 ms and t = 80 ms) , but the two stable 





150 - 


[m/s] 


100 - 













-50 




-1( 



t= 10ms 
t= 50ms 
t= 80ms 
t= 90ms 
t= 100ms 




5 10 
E r [kV/m] 



15 



20 



FIG. 3: (Color online) Ambipolar equation at p = 0.7 for 
several representative times. 




FIG. 4: (Color online) Radial electric field profile for 
representative times. Lines are separated by 1 ms and 
thick lines by 5ms. The starting (t<t cr ) E r (p) is 
positive. 



roots are separated by an unstable root, so E r stays posi- 
tive. Only for when the electron root disappears 
(t = 90ms and t — 100 ms), the system jumps to negative 
E r in tens of ps (the typical time of evolution towards 
ambipolarity, see below). Fig. 4 shows, for the first time, 
the formation and precise evolution of the shear layer 
in TJ-II plasmas. It starts to develop at p = 0.85, ap- 
proximately where the gradient is maximum [9] and then 
propagates inwards and outwards: a speed of the order 
of 1 m/s may be extracted, as measured in Ref. [14]. Let 
us recall that, since we neglect the shear-flow viscosity, 
our simulations are local, hence the evolution equation 
for E r is solved independently at each radial position, 
and only an indirect coupling exists through the n and 
T e gradients. Thus, the speed of the layer propagation 
is determined basically by the evolution time of the local 
collisionality (as we have discussed, n and T e ). Shearing 
rates of the order of 10 5 s _1 can be inferred from Fig. 4, 
so this neoclassical shear might be large enough for play- 
ing a role in the reduction of turbulent transport (the 
linear growth rate of the resistive interchange instability, 
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FIG. 5: (Color online) Density dependence of the 
neoclassical poloidal viscosity during the transition. n cr 
is locally defined as the n at which E r passes through 
(inset: frequency dependence of the neoclassical 
damping for several values of fi p ). 



that is thought to be the dominant one at the edge of 
TJ-II, is also of the order of 10 5 s" 1 [15]. 

The main features of the low density transition are well 
captured by neoclassical transport and the correct pre- 
diction of the shear layer supports the initial ordering 
assumption. If the shear-flow viscosity had been calcu- 
lated, the shear of E r would be somewhat smaller [26], 
but this agreement with the experiment shows that no 
qualitative differences would have been found. 

Generally, this picture is only slightly modified by con- 
sidering the measured turbulent momentum transport. 
Indeed, local measurements of Reynolds stress [6] show 
average events r^s/n^O.l m/s (which would not be per- 
ceptible in the evolution of E ri see Fig. 3) and infrequent 
extreme events of Tjis/n< 10 m/s (which would not fun- 
damentally modify the results). This is however not true 
when the transition is approached from below as can be 
seen in Fig. 3. In this situation (t = 80 ms), the non 
ambipolar neoclassical fluxes display a weak dependence 
on E r around the ambipolar value and large E r excur- 
sions may be caused by turbulent momentum fluxes or 
external forcing (biasing) as is observed experimentally. 

To make the argument more precise we define a neo- 
classical poloidal viscosity as the linear coefficient of 
the difference between the electron and ion radial fluxes 
(shown in Fig. 3) expanded around the ambipolar elec- 
trie field, [T e -Ti}(E r )= -^(E r - E° r ) +0((E r - E° r f). 
Eq. (3) then yields: 



dE r 
~df' 



47r 2 mn(ee • e#) 
-v p (E r -E°)^ 



fi p (E r — E®) 



• Vr) 



(5) 



For the sake of argument we have neglected the (assumed 
slow) variations in n and T^. The coefficient \i v can be 
calculated directly from the data of Fig. 2, and we have 



absorbed n, m, e, ?//, (e# • eg), and constants into v v and 
j r . This scalar viscosity is a combination of the elements 
in the viscosity matrix defined in Ref. [25] (in particular, 
in the limit where the parallel currents are low and the 
temperature gradients small, it is proportional to the L\ 
coefficient). The dependence of v v on n during the transi- 
tion is shown in Fig. 5: it is smaller before the transition 
than after it (as predicted, with a simplified formulation, 
in Ref. [7]) and, more importantly, goes to zero when 
approaching n cr from n < n cr . We now show that the 
behaviour of the neoclassical viscosity provides a simple, 
unified explanation of the observed phenomena (ampli- 
tude of low-frequency, LRCed potential fluctuations, po- 
tential relaxation time and shear-flow susceptibility) that 
accompany the transition. To our knowledge, this is the 
first time that the vanishing of the neoclassical viscosity 
is explicitly shown to be enough to cause the phenomena 
at the transition. 

The characteristic relaxation time in biasing turn-off 
experiments {y~ x in Eq. 5) shows a peak around the crit- 
ical density and decreases for larger density plasmas in 
the ion root [10]. This is reproduced by the curve shown 
in Fig. 5. Similarly, when a low frequency external bias- 
ing is applied, the response electric field is in phase with 
the biasing and its amplitude increases close to the criti- 
cal density [17]. This is to be expected from Eq. 5, for in 
that case E r (t) w E® + v~ l j r (t). Finally, to better dis- 
cuss the observations of LRCs close to the transition [8] 
we Fourier transform Eq. (5): 



iuE r (uj) = - 
\E r (u)\ 2 



1 |j(w)| 2 = A(lo)\j(uj)\ 2 , (6) 



V 2 + OJ 2 



for time scales faster than that of the density ramp, i.e., 
uj> d t \og(E®) ~ d t log(n) ~ 10 Hz. Eq. 6 shows that the 
amplitude of the fluctuations E r (uj) driven by a given 
broadband turbulent forcing is modulated by the 
neoclassical viscosity, which damps fluctuations of fre- 
quencies lower than u p . Below and above the transition 
(Fig. 5 inset), the fluctuations E r (uo) with uj < 10 kHz 
are neoclassically damped. It is only close below n cr 
that the neoclassical viscosity drops, leaving the low fre- 
quency E r fluctuations (which are expected to display 
higher LRC) undamped. A higher order shear-flow vis- 
cosity [26] is usually included in Eq. (3), whose main ef- 
fect is to smooth the radial variations of the E r . For the 
present discussion, this term might only be important in 
the inner side of the shear layer, where the leading order 
viscosity goes to zero, and will set a low (though not null) 
effective viscosity and damping rate. The experimental 
measurement [10] of a peaking factor of about 5 in v~ x 
gives an upper limit to these second order effects. 

Note that the vanishing of the viscosity is a conse- 
quence of the transition from electron to ion root and 
the above phenomena are not to be expected in the re- 
verse transition [12, 16]. 
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In conclusion, we have described, by solving the drift 
kinetic equation, the formation of the shear layer in TJ- 
II. Even though such a first-principles theoretical calcula- 
tion is by itself an important result, we have additionally 
advanced in the understanding of relevant physical phe- 
nomena associated to the transition, and also observed 
in other devices. They are essentially related to the am- 
plification of dynamical electric fields as a result of the 
reduction or vanishing of the neoclassical viscosity. 
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